source("Rep RDD.R", print.eval=TRUE)
### saving p-values ###
IT1a <- rdlamp1$p[1]
IT1b <- rdlamp2$p[1]
IT2a <- rd13IT$p[1]
IT2b <- rd13bIT$p[1]
IT3a <- rd15IT$p[1]
IT3b <- rd15bIT$p[1]

cc1a <- rd1$p[1]
cc1b <- rd1b$p[1]
cc2a <- rd2$p[1]
cc2b <- rd2b$p[1]
cc3a <- rd3$p[1]
cc3b <- rd3b$p[1]
cc4a <- rd4$p[1]
cc4b <- rd4b$p[1]
cc5a <- rd5$p[1]
cc5b <- rd5b$p[1]
cc6a <- rd6$p[1]
cc6b <- rd6b$p[1]
cc7a <- rd7$p[1]
cc7b <- rd7b$p[1]
cc8a <- rd8$p[1]
cc8b <- rd8b$p[1]
cc9a <- rd9$p[1]
cc9b <- rd9b$p[1]
cc10a <- rd10$p[1]
cc10b <- rd10b$p[1]
cc11a <- rd11$p[1]
cc11b <- rd11b$p[1]
cc12a <- rd12$p[1]
cc12b <- rd12b$p[1]
cc13a <- rd13$p[1]
cc13b <- rd13b$p[1]
cc14a <- rd14$p[1]
cc14b <- rd14b$p[1]
cc15a <- rd15$p[1]
cc15b <- rd15b$p[1]




### new dataset of p values ###
p <- c(IT1a,IT1b,IT2a,IT2b,IT3a,IT3b,cc1a,cc1b,cc2a,cc2b,
          cc3a,cc3b,cc4a,cc4b,cc5a,cc5b,cc6a,cc6b,cc7a,cc7b,cc8a,cc8b,cc9a,
          cc9b,cc10a,cc10b,cc11a,cc11b,cc12a,cc12b,cc13a,cc13b,cc14a,cc14b,
          cc15a,cc15b)

p <- as.data.frame(p)




### plot p values in comparison to uniform distribution ### 
theme_base2 <- 
  theme_minimal(base_size=20)  + 
  theme(legend.position=c(0.15, .9),  legend.key.size = unit(1,"line"),
        axis.text=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),
        plot.title = element_text(size=20, hjust= 0.5))


ggplot(p, aes(sample = p)) + 
  stat_qq(distribution = stats::qunif) +  
  geom_abline(linetype = "dashed", size = 1.2) + 
  geom_point(stat="qq", distribution = stats::qunif, size= 5) +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "red3", size = 1.2) + 
  geom_vline(xintercept = 0.1, linetype = "dashed", color = "blue", size = 1.2) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "red3", size = 1.2) + 
  geom_hline(yintercept = 0.1, linetype = "dashed", color = "blue", size = 1.2) +   
  theme_light(base_size = 35) + ylab("p-values") + 
  xlab("Uniform Distribution") +
  ylim(0,1) +
  xlim(0,1)
